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1 INTRODUCTION 



' The r-mode i nstab ility iAnderssorJ I199S 

. iFriedman fc MorsinkI Il998il in rotating neutron stars 
' has attracted much attention due to its possible role in 
, limiting the spin rates of neutron stars and the possibility 
I of producing gravitational waves detectable by LIGO II. 
However, the significance of this instability depends strongly 
on the maximum amplitude to which the r-mode can grow. 
A number of damping mechanisms of the r-mode have 
been examined, but no defi nite conclusion of the maximum 
amplitude can be made (seelAndersson fc Kokkota^ jgOO^): 
IFriedma n fc LockitchI ^2001^ for reviews: see also lAnderasonl 
(^003) and references therein for a discussion of recent 
development s) . 

In Ref. llGressman et al.ll2002^ ■ we reported that large 
amplitude r-modes are subjected to a hydrodynamic insta- 
bility that worked against the gravitational-radiation driven 
instability. Our numerical simulations showed that, for an 
r-mode with a large enough amplitude (e.g., with a dimen- 
sionless amplitude parameter a of order unity), the hydro- 
dynamic instability operated in a timescale much shorter 
than that of the gravitational-radiation driven instability. 
This hydrodynamic instability when fully developed will 
cause the r-mode to decay catastrophically. As a result, the 
neutron star changes rapidly from the uniformly rotating 
r-mode pattern to a complicated differential rotating con- 
figuration. However, we have not been able to identify the 
hydrodynamic mechanism responsible for this catastrophic 
decay. 

In this paper we establish that this catastrophic decay is 
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due mainly to a particular three-mode coupling between the 
r-mode and a pair of fiuid modes. We identify the daughter 
modes and show that their frequencies satisfy a resonance 
condition: the sum of their frequencies is equal to the r- 



mode frequency. The importance of three-mode coupling in 
the r-mode evolution was first pointed out bv Arra s et alJ 
(2003). By comparing the fully nonlinear simulation with 
the second-order perturbation equations, we determined the 
responsible mode coupling strength k. The value of k im- 
plies that this particular three-mode coupling alone gives a 
rather stringent limit to the usual picture of a neutron star 
with a fluid flow pattern dominated by the r-mode: Such a 
picture can only be valid with the dimensionless r-mode am- 
plitude a smaller than or of the order 10~^. Beyond that, 
the r-mode must be considered to be strongly coupled to 
other fiuid modes. Instead of the gravitational radiation re- 
action, it is the fiuid coupling that dominates the dynamics 
of the r-mode when its amplitude grows to 10"^. Its evolu- 
tion is strongly affected by the fluid coupling in a timescale 
much shorter than its growth time due to the gravitational- 
radiation driven instability. In particular, at that point the 
steady rise in amplitude of the r-mode due to gravitational 
radiation reaction will turn into a catastrophic decay. This 
is the main conclusion of the paper. 

What happens after the catastrophic decay? As the 
post-decay fiuid fiow involves couplings between a large 
number of short wavelength modes beyond the resolution of 
our simulation, we cannot rule out the possibility that the 
r-mode can regain energy from these fiuid modes and be re- 
vived within a hydrodynamic timescale (for a detailed study 
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of such "bouncing" phenomena see lBrink et alJ l|200^) be- 
fore the star returns to uniform rotation. However, we do not 
see any hint of such bouncing in our simulations with an r- 
mode amplitude of order unity. This observation together 
with the consideration that a large number of fluid modes 
are involved in the catastrophic decay, lead us to conjecture 
that there is no such "bouncing" of the r-mode amplitude 
after the catastrophic decay. 

We further note that there could be fluid modes with 
wavelengths shorter than what our numerical simulation can 
resolve but with couplings to the r-mode strong enough that 
pl ace a stronger limi t on the r-mode amplit ude, as argued 
bv lArras et"all (l2003f) and lBrink et al.l Jiooi). However, the 
bound we establish in this paper already has significant im- 
plications on the astrophysical importance of the r-mode 
instability. 

The study in this paper makes no particular assump- 
tion on the neutron star model. The star is described sim- 
ply as a self-gravitating, perfect-fluid body with a generic 
equation of state (a polytropic EOS); no dissipative mech- 
anism is assumed. We have seen the same hydrodynamic 
phenomenon in simulations with stellar models of different 
central densities and polytropic indices. This suggests that 
the hydrodynamic mechanism being studied is quite generic. 



2 NUMERICAL RESULTS 
2.1 Nonlinear mode coupling 

We solve the Newtonian hydrodynamics equations for 
a non-viscous, self-gravitating fluid body in the pres- 
ence of a current quadrupole post-Newtonian radiation 
reaction ('Gressman et al.l l2002f) . We refer the reader to 
[Grossman ct al. (2002) for the neutron star model used in 
this study; and a detailed discussion on how we set up a 
large amplitude r-mode as an initial state for simulation. To 
set the stage for the study in this paper, we show in Fig. 
0the evolution of the r-mode amplitude a vs. time for the 
case where the initial nonlinear r-mode configuration has 
amplitude a = 1.6. We compare the results for three differ- 
ent grid resolutions 129^, 161^, and 193^. The inset in Fig. 
Q shows the details of the evolution in the time interval (15 
ms-25 ms). Fig. Q shows that a starts off slowly decaying 
until t « 40 ms (in the 193^ simulation) at which point it 
decays catastrophically. 

To investigate what hydrodynamic effect is responsible 
for the catastrophic decay of a, we divide the slowly decaying 
phase into a few time intervals, and study the Fourier spec- 
tra of the velocity fields. In Fig. |21(a), we plot the Fourier 
transform of the axial velocity along the x axis at a; = 2.56 
km inside the star for the 129^ result plotted in Fig.0 The 
dashed line is the Fourier spectrum of the earlier part of the 
evolution (13 ms-35 ms), while the solid line represents that 
of the later part (35 ms-50 ms). We see that the dashed line 
has only one large peak at the r-mode frequency (ro = 0.93 
kHz). This is compared to the spectrum at the later time slot 
showing many peaks at different frequencies. In particular, 
we see a significant peak at 0.53 kHz. 

Fig. 121(b) shows the corresponding Fourier transforms of 
the rotational velocity along the x axis at the same point 
as in Fig.|5|(a). In contrast to Fig.|5|(a), there is no peak at 
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Figure 1. Evolution of a with initial value 1.6 with three dif- 
ferent resolutions. The inset shows the details of the evolution 
between 15 ms-25 ms. In the 193^ simulation, the neutron star is 
represented by a total number of fs! 4.4 X 10^ grid points. 

the r-mode frequency ro in the initial time interval (dashed 
line). This shows that during the initial evolution the large 
amplitude r-mode maintains its character as given in the 
linear theory: no component on the equatorial plane. In 
the later time interval (solid line), we see various peaks, 
and in particular a significant one at 0.40 kHz. It is also 
interesting to note the appearance of the peak at the r-mode 
frequency ro in the later time evolution. This is consistent 
with the formation of the large-scale vortex of the r-mode 
flow pattern du r ing th e late time evolution as reported in 
iGressman et alj i2002l) . The r-mode flow pattern changes 
rapidly to a vortical motion on the equatorial plane at the 
time when a collapses rapidly. 

Figs. |5| (a) and (b) demonstrate that the r-mode can 
leak energy to other fluid modes through nonlinear hydro- 
dynamic couplings, and in particular leaks energy to the 
fluid modes at frequencies 0.53 (henceforth the "^4" mode) 
and 0.40 kHz (the "B" mode). Note that the sum of the 
frequencies of the A and B modes exactly equals to the r- 
mode frequency (ro — 0.93 kHz). This relation suggests that 
the r-mode leaks energy to this pair of modes via a three- 
mode couplingrnechanism satisfying the resonance condi- 
tion liCraiklll985h 

f.A + fB = fR, (1) 

where the f's are the mode frequencies. 

How do the coupled fluid modes grow? Instead of just 
two time slots, we can further divide the slowly decaying 
phase of a into several subintervals and perform Fourier 
transform on each of them, and monitor the Fourier ampli- 
tudes of different time intervals. In Fig.|3 we plot the Fourier 
amplitude of the A mode vs time for the case q = 1.2 with 
resolution 193^. In the figure, the circles connected by a solid 
line are numerical data. The first data point corresponds to 
the Fourier amplitude obtained in the time interval 10 ms- 
25 ms. The time value of the data point (17.5 ms) is taken 
to be the middle of this interval. The second data point cor- 
responds to 15 ms-30 ms, etc. The dashed line is fitted to 
the exponential function / — Aoexp{{t— 17.5)/9.35), where 
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Figure 2. (a): Fourier spectra of along the x axis (at x = 2.56 km) in the evolution of Fig. 1 at two time slots in the 129'^ simulation, 
(b): Fourier spectra of for the same case. 



t is in millisecond and Aq is the Fourier amplitude obtained 
in the interval 10 ms-25 ms. Fig. |3 shows that the A mode 
grows exponentially with a time constant « 9.35 ms during 
the slowly decaying phase of a. The B mode grows expo- 
nentially with a time constant about 7 percents larger as 
determined numerically. 

Our results thus show that the coupled fluid modes are 
unstable and grow exponentially due to their couplings to 
the r-mode. An important point that we observed in all cases 
studied is that when the amplitudes of these two coupled 
daughter modes grow to a point comparable to that of the 
r-mode, the catastrophic decay of the r-mode sets in. 



2.2 Determination of the Decay Rate of a 

In the slowly decaying phase of a, there are two mechanisms 
that lead to the decay of a in our numerical simulation: 

(i) numerical damping arising from flnite differencing error, 
which depends on the resolution used in the simulation; and 

(ii) physical damping due to the nonlinear interactions of 
the r-mode to other fluid modes, which depends on a at 
that moment. Fig. shows that numerical damping is neg- 
ligible during the early part (15 ms-30 ms) of the evolution 
at high enough resolutions. The 161^ and 193^ results agree 
to high accuracy. We can thus regard the rate of change of 
the r-mode amplitude da/dt in the 193^ simulation as a rea- 
sonable approximation to its physical value (when Aa; 0) 
due to the particular three-mode coupling seen in our simu- 
lation. A linear curve-fit yields 

da 
Hi 



— = -2.9 X 10"'^ (ms)" 



(2) 



2.3 Determination of the waveforms 

To further characterize the daughter modes, we would like 
to determine their waveforms in addition to their mode fre- 
quencies. However, this turns out to be non-trivial. 
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Figure 3. Exponential growth of the Fourier amplitude of the A 
mode for the case a = 1.2 with resolution 193^. 



In Fig. 131(a), 



show the velocity profile for the 



129"^ simulation shown in Fig. along the x axis aX t — 
28.88 ms. We note that this is basically the velocity profile 
of the m = 2 r-mode. To extract the waveform of the 
A mode, we need to subtract out the r-mode waveform. 
However, the r-mode waveform is known only in the linear 
approximation, while the r-mode being studied in this paper 
is not small and there can be non-linear correction to the r- 
mode waveform. To do the subtraction, we note that the 
non-linear correction should not change the symmetry of 
the r-mode, and in particular that its is symmetric with 
respect to reflection about the origin on the x axis. We thus 
deflne the quantity Aw^ := {v^{x) — {—x))/2 on the x 
axis. In Fig. 21 (b), we plot Ad^ in the core region of the 
star at three different times for the 129^ simulation shown 
in Fig. [H t = 28.88 ms, t = 29.82, and t = 30.76 ms. The 
latter two time values differ from the flrst one respectively 
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Figure 4. (a): The velocity profile along the x axis at t = 28.88 ms for the 129"^ simulation shown in Fig.0 (b): Waveforms of the 
A mode for the same case. Plotted are the profiles of Ai)^ (see text) along the x axis at three different times. The latter two time values 
differ from the first one respectively by half and one oscillation period of the A mode. 



by half and one oscillation period of the A mode (recall 
that /a = 0.53 kHz). We see explicitly that the waveform 
oscillates with the frequency of the A mode. The two lines 
aX t = 28.88 and 30.76 ms match pretty well, except in the 
region near the surface of the star where the r-mode velocity 
is large. This suggests that the extracted waveform is indeed 
the eigenfunction of the A mode. Note that the amplitude 
of the A mode is less than 1% of that of the r-mode in this 
early slowly-decaying phase. The extraction would not be 
possible without taking advantage of the symmetry. 

Next we turn to the B mode, which is a mode with fluid 
flow in v'" . One might expect the waveform extraction to be 
easier, as the r-mode has no v"" component. To begin, we 
subtract the background rotation profile from the velocity 
component n'' by defining An^ ~ v^{x) — Vq{x) along the 
X axis, where is the profile at the starting point of the 
slowly decaying phase of a. In Fig.|5] we plot Av^ along the 
X axis at various times in the slowly decaying phase for the 
same 129^ simulation as before. We see a small differential 
rotation pattern growing up to an amplitude a few percent 
of that of the r-mode. Four features of this fluid pattern 
are: 1. This differential rotation mode has a long period, 
much longer than that of the B mode. 2. The amplitude of 
this mode saturates at around 35 ms. 3. The energy of this 
mode is small comparing to the r-mode. 4. This mode is 
antisymmetric with respect to the origin along the x axis. 

We uncovered a differential rotation pattern unrelated 
to the A and B modes. What is its significance? We make 
two observations: 1. One might suspect that this differential 
rotation pattern is part of the non-linear r-mode profile. It 
has been argued that differential ro tation can be induce d 
by second -order r-mode motion (e.g.. lRezzolla et all i200Cl) : 
ISa l|2003)) in perturbation analyses. However, the differen- 
tial rotation_2attemjhow^ 1^ does not match those 
of ('RczzoUa e t"alll2000t Isill2004l) . Further, we see that the 
differential rotation given in Fig. I^is growing in amplitude 
during the slowly decaying phase of the r-mode. This sug- 
gests that it may not be part of the r-mode. 2. This differen- 
tial rotation pattern is different from the one that developed 
during the catastrophic decay of the r-mode as we reported 



previously iCressman et alJl2002l) . Further investigations of 
this unknown mode {U mode) are needed. 

To extract the B mode, we subtract the antisymmetric 
pattern of this "[/ mode" and find a symmetric (with respect 
to reflection about the origin) component with a period 
the same as that of the B mode. It is smaller than the U 
mode at earlier time, but overtakes the U mode in amplitude 
at around 40 ms. It keeps on growing exponentially to an 
amplitude comparable to that of the r-mode at the point of 
the catastrophic decay. In Fig. |S|we show the waveforms of 
the B mode (obtained by subtracting the U mode pattern 
from the profile of Av^) at around the time it becomes the 
dominant pattern. The two time slices 40.02 ms and 41.22 
ms are separated by half of the oscillation period of the B 
mode, and indeed we see that the two lines are off by vr in 
phase. 

To summarize, we found an antisymmetric waveform in 
for the A mode, and a symmetric waveform in for the 
B mode. The wavelength of both modes are about 4 km in 
the central region of the star. Next we want to investigate if 
the symmetry properties of the waveforms of the modes are 
expected. We note that under a 7r-rotation (i.e., <j!> — > <^ + tt) 
on the equatorial plane, the Cartesian components of the 
velocity field of a toroidal fluid mode^ transform according 
to t;" ^ (-l)™+lu^ and ^ i-l)"'v\ 

where m is the azimuthal mode number (see Appendix IXll. 
This implies that the profiles of and along the x axis for 
an m = — 1 mode are symmetric with respect to refiection 
about the origin; while the profile of for an m = — 1 mode 
is antisymmetric. We note that the A and B modes having 



^ Toroidal modes are characterized by fluid displacement vectors 
proportional to r X Vyjm, where Yj^ is the standard spherical 
harmonics. The r-mode is known to be a torodial mode. We as- 
sume that the two daughter modes seen in our simulations are also 
torodial modes (to leading order) because of their low frequencies 
compared to spheroidal modes (e.g., the /- and p-modes). The 
low frequency (spheroidal) g-modes do not exist in our isentropic 
simulations. But note that the inertial modes of isentropic stars 
in gen eral have a "hybrid" nature (see .Andersson Kokkota3 
J2OO J) for discussion). 
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Figure 5. DifTcrential rotation pattern Av^ (see text) for the 
same simulation as Fig. Bl 
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Figure 6. Waveforms of the B mode (see text) along the x axis 
at two different times (with a time interval half the oscillation 
period of the B mode) for the same simulation as Fig. |1] 



m = —1 is consistent with the picture of three- mode cou- 
phng. The r-mode has m = 2 and the conservation of angu- 
lar momentum requires that the azimuthal mode number m 
of the three couple d modes must satisfy mu + -|- ms = 
l|Arras et aljl2nn,t ). 



3 SECOND-ORDER PERTURBATION 
ANALYSIS 

3.1 Three-mode coupling 

The numerical results indicate that the dominant hydrody- 
namic contribution to the decay of the r-mode is due to its 
coupling to the A and B modes. To gain further insight into 
the hydrodynamic evolution during the early part of the nu- 
merical evolution, we compare the fully nonlinear simulation 
to a second-order perturbation analysis. In particular, we 
want to estimate the coupling strength of the three-mode 
coupling found, making use of the decay rate of a deter- 
mined in Eq. The second-order perturbation analysis 
enables us to extract information about the mode coupling 
independent of the specific value of a used in the numerical 
simulation. Through this we will then be able to determine 
the behavior of the r-mode at a values smaller than that 



used in the simulations. Note that the following analysis 
only makes use of information obtained in the slowly decay- 
ing phase of a, where a second-order perturbation analysis 
is assumed to be valid. In particular, it does NOT depend 
on the existence of a catastrophic decay. The perturbation 
description is NOT expected to be valid when the evolution 
is getting near to the point of the catastrophic decay and 
beyond. 

In a reference frame co-rotating with a rotating star, the 
second-order perturbation equations fo r a system of three 
modes in resonance can be written as dSchenk et alj|200ll : 
lArras et a^]|2003^ 



<1R , ■ 
at m 



dqA 

~dr 

dqs 
dt 



+ iujAqA = -JAqA + iuJAHRABqnqB, 



+ ilJJBC 



-iBqB +ii^Bl^RAB<lRqA, 



(3) 



(4) 



(5) 



where the g's represent the dimensionless complex ampli- 
tudes of the modes; the oj's are the rotating frame angu- 
lar frequencies of the modes; tr is the growth time of the 
r-mode due to radiation reaction; krab is the coupling 
strength; Z* denotes the complex conjugate of Z; ^a and 
7B are the damping coefficients of the two coupled modes 
representing the outflow of energy from the three mode sys- 
tem. 

The system of equations lIHJ-lO gives rise to a resonant 
coupling for modes satisfying ua + uib + i^R ~ 0. Note that 
this is not inconsistent to the resonance relation Q satisfied 
by the (inertial frame) frequencies of the coupled modes. 
The rotating frame angular frequency of a fluid mode u is 
related to the inertial frame angular frequency uio hy uj — 
uiQ + mQ, where il is the angular velocity of the star. For 
the star model (n = 4.87 (ms)~^) and m = 2 r-mode with 
the inertial frame frequency ro = 0.93 kHz being studied, 
we have ujr = 3.9 (ms)~^. Note that we deflne loq to be 
negative for the r-mode (see iGressman et~a l. (2002)). The 
rotating frame angular frequencies of the m = — 1 daughter 
modes are respectively ua = —1.54 (ms)~^ ~ —OAlor and 
ujb = —2.36 (ms)~^ ~ —O.Qujr, suggesting that the three 
modes indeed satisfy the resonance condition Suj = u>a + 
+ (^fl ~ 0. 

The me aning of the 7's us ed h ere is somewh at different 
from that of lArras et ail J2003h and lBrink et al.1 (12005.) . The 
neutron star EOS used in this paper is that of a perfect 
fluid with no viscosity. The 7's describe the couplings of 
this particular three-mode system to the other fluid modes 
of the neutron star and is a function of the amplitude of these 
other modes. The description of Eqs. (I21l-(|Hll, is meaningful 
only when the 7's are not changing in a time scale short 
compare to that of the dynamical time scales of the A, B 
and r modes, a point we will return to later i n the paper. 

A general expression for krab is given in lSchenk et alJ 

i200ll) . In writing the coupled equations above, we omit the 
coupling terms proportional to kraa and krbb, which ar e 
zero due to a z-parity selection rule dSchenk et alJlioOll) . 
This selection rule states that habc = for an odd number 
of 2-parity odd modes. Since the r-mode has odd z-parity, 
i^raa = regardless of the z-parity of mode A. In order that 
Krab 7^ 0, modes A and B should have different z-parity. 
We see in our simulations that one of the coupled modes (A) 
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appears only in the Fourier spectrum of , while the other 
one [B) appears only in the spectrum. 

No te that th e definition of krab used by us is the 
same as lArras e t al. ( 2003 ), b ut different from tha t used by 
iBrink et al.1 ll2005l) fand also ISchenk et alJ ll200ll) '). In our 
notation, krab is the ratio of the nonlinear interaction en- 
ergy to the energy of the mode (both at unit amplitude) . 
The c onversion from our k to the one used bv lBrink et alJ 
^2005^ is K ^ EunitK, where E^nit = \MB?Q^ is the mode 
energy at unit amplitude in the normalization used by them. 

In order to relate qn to the nonlinear r-mode amplitude 
defined in our numerical simulations we take [(/^(t)! — ctit) 
(the critical amplitude ctcrit given below applies as long as 

is proportional to a). Using Eq. (^ and its complex 
conjugate, it can be shown that a satisfies 



da a ujR 

= \ lm(Ki{As(7A(7s(7fl), 

TR a 



dt 



(6) 



where \m[Z) denotes the imaginary part of Z. In the slowly 
decaying phase of a, \qA\ — |gsl = g is less than 
(with the catastrophic decay sets in at g ~ a). We obtain 
the following order-of-magnitude estimate for the coupling 
strength \krab\: 



1 



\hr.ab\ 



LOrQ^ 



da 
It 



a 

TR 



3 X 10" 



(7) 



where we have used a = 1.6, da/dt — —2.9 x 10^ (ms)^ , 
ijR = 3.9 (ms)~^, TR = 40 s, and q « a. Note that using a 
smaller value of q will only increase the value of |khas|. In 
this sense, the \hrab\ given by Eq. JTJ is a lower bound. A 
larger |Kfl,4s| gives a tighter bound to the r-mode amplitude 
at the point of catastrophic decay Ocrit as shown below. 

To check that the value of \krab\ does not depend 
sensitively on the initial value of a used in the simula- 
tion, we have repeated the above calculation with a dif- 
ferent initial r-mode amplitude a — 1.2. In that case, we 
obtain da/dt — —1.7 x 10""^ (ms)~^ during the slowly de- 
caying phase of a. The value of |KijA_B| changes only by a 
few percents (within numerical error) comparing to the case 
a = 1.6. 

To complete the description of the system using Eqs. 
©-lIHJ, one would also like to determine from the numerical 
simulation the values of "/a and 7_b. In our present setup, 
7a and 73 determine the rate of energy leaking out of the 
three-mode system through coupling to other fluid modes. 
However, due to the constraint of computational resources 
we were not be able to obtain the precise values of 7a and 
7b since they correspond to a timescale orders of magni- 
tude longer than 10 ms (the timescale of the growth of the 
daughter modes in our simulations) , in the regime where the 
daughter modes A and B are small compared to the r-mode. 

In our simulations, when the amplitudes of the two 
daughter modes become comparable to that of the r-mode, 
the perturbation equations (|^-l|^ are no longer valid. The 
evolution becomes fully nonlinear, and the r-mode decays 
catastrophically at that point. 



3.2 Amplitude of the 7?-mode at the catastrophic 
decay 

Now we consider how an r-mode would evolve from a small 
initial amplitude to the decay point in a star described by 



a perfect-fluid EOS, based on the fully nonlinear numerical 
simulation and information from the second-order pertur- 
bation analysis. We assume that during the early evolution 
the r-mode couples mainly only to the two resonant modes 
{A and B). In this phase, the system can be described by 
Eqs. (|2J-(0. The A and B modes are negligible initially; the 
r-mode grows due to radiation reaction. When a becomes 
large enough, \qA\ and grow rapidly. The point at which 
the two daughter modes catch up in amplitude with the r- 
mode can be estimated by setting da/dt — in Eq. Q. 
In particular, the amplitude of the r-mode at that point is 
given by 

(8) 



Qcrit 



0.02, 



(^rtr\krab 



where we have taken \qA\ ~ \qB\ ~ Ocrit. 

At this point, the usual picture of a neutron star with 
fluid flow dominated by an r-mode pattern is no longer valid. 
It is somewhat surprising that a single three-mode coupling 
of purely hydrodynamic nature could put a limit as strong 
as a ~ 0.02. To the best of our knowledge this is the first ex- 
plicit demonstration that an r-mode cannot be the dominant 
feature of a rapidly rotating neutron star with an amplitude 
beyond a ~ 0.02. 

The critical amplitude Qcrit in Eq. Q is defined to be 
the amplitude of the r-mode when nonlinear hydrodynamic 
couplings dominate the gravitational radiation reaction in 
driving the r-mode evolution; it is the time when the daugh- 
ter modes catch up in amplitudes with the r-mode. It should 
be noted that this critical amplitude is different from the so- 
called par ametric instability threshold (see, e.g., Eq. (2) of 
iBrink e~al. (2005) ), which is defined to be the amplitude of 
the parent mode at the point when the two daughters be- 
come unstable due to nonlinear coupling and start to grow 
exponentially. The scalings of the two amplitudes are differ- 
ent. In particular, the parametric instability threshold does 
not depend on the growth time tr of the r-mode, while our 
critical amplitude Ocrit is inversely proportional to tr. 



3.3 Beyond the Critical Point 

What happens after the point a = Qcrit is reached? Our nu- 
merical simulation shows that when the two daughter modes 
catch up in amplitude with the r-mode near Qcrit, the follow- 
ing happen: 1. Many other fast growing fiuid modes enter the 
picture. 2. The fiuid fiow can no longer be approximately de- 
scribed by the three-mode perturbation system Eqs. 
3. The r-mode decays catastrophically. 4. After the catas- 
trophic decay, the neutron star has a complicated differential 
rotating fiow pattern with short length scales, and 5. there 
is no sign of the r-mode re-gaining its energy from the many 
fiuid modes it directly or indirectly coupled to. 

One important implication of this is: The maximum am- 
plitude Qmax an r-mode can get to through the gravitational- 
radiation driven instability is given by it. How- 
ever, there is a caution: our simulation is carried out with a 
of order unity but not directly at q ~ 0.02. There is a possi- 
bility that even though the fiuid fiow is no longer dominated 
by that of the r-mode after the two daughter modes catch 
up at Qcrit ~ 0.02, the three-mode perturbation system Eqs. 
(Ol-© could remain a valid description of the neutron star. 
The perturbation system can fail to represent the fiow of the 
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star for two reasons: (i) the mode amplitudes become large, 
so that higher order perturbation terms in the equations can- 
not be neglected, or (ii) the effects of modes other than the 
three modes studied become significant (this shows up also 
as the 7's are changing on a short time scale); such a situa- 
tion is strongly suggested by our simulation. Unfortunately 
simulations of timescales long enough to accurately capture 
the behavior of r-modes with a small a value require a com- 
putational resource beyond what is available to us, and we 
cannot investigate the catastrophic decay at small a value 
in full. 

In case the second-order perturbation equations 
are still valid at the point a = Oait, what could happen sub- 
sequent to the catastrophic decay? In particular, can the am- 
plitude of the r-mode evolve to a value substantially larger 
than ctcrit? The subsequent evolution described by Eqs. Q- 
JSJ depends strongly on the relative values of 7^, 7_g, th, and 
th e frequency detuni ng 5uj = uia + <jJs + i^R- As discussed 
in lArras et alj i2003fl . a stable equilibrium solution of Eqs. 
0-l|5Jl exists if both the damping and detuning are large: 
7A + 7s > r^^; \3uj\ > {■ja + 7b)/2. Otherwise, the system 
is unstable. (Note that in our study 7a and 7s represent 
the rate of energy leaking out of the three-mode system to 
other fluid modes of the neutron star, which is described by 
a perfect fluid EOS with no viscosity.) 

In the case where a stable equilibrium solution exists, 
the system will settle into a steady state with the mode am- 
plitudes being roughly the same. The maximum amplitude 
of the r-mode Umax would be given by Ofcrit- On the other 
hand, in the unstable case Eqs. (|^-(|KJ predict that after the 
amplitudes of the three coupled modes become comparable, 
the amplitudes of the three modes will fluctuate in a short 
timescale with rapid energy transfer between them, and the 
average amplitude of all three modes could grow until the 
second-order perturbation analysis becomes invalid. How- 
ever, such oscillation of energy between the three modes is 
possible only if their phase relations are well preserved in the 
evolution, a condition we believe may not happen in realis- 
tic nonlinear hydrodynamic evolutions: at a « Qcrit a large 
number of other fluid modes enter into the picture and affect 
the fluid evolution significantly. The energy originally in the 
r-mode is transferred to many fluid modes. The involvement 
of a large number of modes suggests that the oscillation of 
energy back to the r-mode in a short dynamical timescale 
is unlikely. As an example, in a chain of A*' coupled non- 
linear oscillators an equipartition of energy among the cou- 
pled modes is achieved if the energy of the initial large scale 
mode is larger than a certain thre shold, which tends to zero 
sufficiently fast at increasing N dCasetti et al.lll99'iD . This 
suggests that, after cascading its energy down to a sea of 
coupled modes, a generic initial large scale mode cannot re- 
gain its energy within dynamical timescale. We believe that 
this is also the case for nonlinear couplings of the r-mode to 
other fluid modes in rotating neutron stars. Hence, we con- 
jecture that such oscillation of energy cannot occur beyond 
the catastrophic decay and the amplitude of the r-mode is 
bounded by Qcrit. 

Nevertheless, we do not rule out the possibility that the 
global r-mode could grow again in a secular timescale (~ 40 
s) due to the driving of gravitational radiation reaction. We 
have extended some of our simulations for a timescale of 10 
ms after the catastrophic decay. However, we do not see any 



evidence of the regrowth of an unstable r-mode, even with an 
artiflcially enlarged radiation reaction force. One possibility 
is that the strong differ ential rotation d eveloped during the 
catastrophic decay feee iGressman et^l. (2002)) may act to 
stabilize the r-mode. If this is the case, it will take even a 
much longer time, until the differential rotation is damped 
out by viscosity (not modelled in our simulations), before the 
r-mode can grow again. To model such regrowth scenario, 
if it occurs at all, would be too expensive computationally 
for present-day technology. This is still an open issue and 
remains to be investigated in the future. 

4 COMPARISON WITH OTHER WORK 

There are two other numerical simulations of nonlinear r- 
modes, both finding that, la rge amplitude r-modes can exist 
for a long period of time. IStergioulas fc Fonli (1200111 per- 
formed relativistic simulations of the r-mode on a fixed 
spacetime background. They started with a large r-mode 
perturbation in a relativistic rotating star and evolved the 
system for about twenty rotation periods. They found no ev- 
idence of mode saturatio n unless the r-mode am plitude was 
much larger than unity. iLindblom et alJ (1200 jl carried out 
numerical simulations of the growth of the r-modes driven 
by the current quadrupole radiation reaction in Newtonian 
hydrodynamics. To achieve a significant growth of the r- 
mode amplitude in a reasonably short computational time, 
they multiplied the radiation reaction force by a factor of 
4500, decreasing the growth time of the r-mode from about 
40 s to 10 ms. The r-mode amplitude a grew to ~ 3.3, be- 
fore shock waves appeared on the surface of the star and the 
r-mode collapsed. They suggested that the nonlinear satu- 
ration amplitude of the r-modes might be set by dissipation 
of energy in the production of shock waves. 

We note that, with the artiflcially larg e radiation reac- 
tion, the results in dLindblom et alj(2001r) assume that no 
hydrodynamic process takes energy away from the r-mode 
in a timescale between 10 ms and 40 s (the artiflcial growth 
time and the actua l physical growth time , respectively). In 
our previous work dCressman et alJl2002^ ■ we showed that 
r-modes of amplitudes unity or above were destroyed by a 
catastrophic decay within this time period. During the pro- 
cess the r-mode motion changed rapidly to a vortical motion 
and strong differential rotation developed. The mechanism 
leading to this catastrophic decay is studied in this paper. 

Another approach to study the nonlinear couplings be- 
tween the r-modes and other fluid modes is via perturbation 
theory. In this analytic approach, one assumes that the sys- 
tem is weakly interactive; and the r-mode couples to other 
fluid modes via the lowest order couplings so that higher or- 
der effects can be ignored. Based on a second-order pertur- 
bation theory developed bv jSchenk et a l. ( 2001), Arras et all 
(2003) proposed a turbulent cascade scenario in which the 
energy of the r-mode leaked to a large number of short wave- 
length inertial modes due to three-mode couplings. The cou- 
pled modes are in turn damped by viscosity. They argued 
that the rotating frame r-mode energy could be limited to 
E, — mode ~ 10~®-Eunit for a millisecond rotating star, where 
£unit ~ ^MR^Q,^ . In their normalization, the r-mode ampli- 
tude (^1 in their notation) was related to the mode energy 
by i5r-modc = 2i5unit^i. Note that the normalization of the 
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such that Al ~ 0.09a (see f 
Hence, in our notation, the maximum r-mode amphtude 
suggested by them is a R:i 8 x 10~^. 

In a recent second-order perturbative work. lBrink et aP 
i2005fl have studied the evolution of the unstable r-mode 
in a nonlinear network coupled with about 5000 inertial 
modes. They considered a uniform density, incompressible 
uniformly rotating star. With this simplified model, they 
calculated the eigenfunctions of the inertial modes and the 
three-mode coupling coefficient k analytically. Their study 
suggested that, while the behavior of the r-mode was com- 
plicated and depended sensitively on the damping effects 
of the coupled modes, the r-mode amplitude would in gen- 
eral be limited to a small amplitude: ~ 10"'* in their 
notation. Note that the normalization convention used by 
[Brink ct al. ( 2005) is such that the r-mode energy is given 
by iSr-modc = -EunitC^. This lu tum gives the relation 
Ca = %/2j4i ~ 0.13a, and thus the r-mode is limited to 
a « 8 X 10~* in our notation. 

In the perturbation studies of iBrink et al.l feOOST) , they 
found that beyond the critical point acrit, energy might os- 
cillate between the r-mode and the fluid modes it coupled 
to. After the catastrophic decay, the r-mode amplitude could 
re-grow to a value comparable to that just before the decay 
in a short timescale. In our numerical simulations with a of 
order unity, we do not see any such oscillations, instead the 
r-mode energy cascades into a large number of fluid modes 
beyond Ocrit. We conjecture that such oscillations are not 
possible in fully nonlinear evolutions and realistic neutron 
stars. This conjecture in turns implies that the r-mode am- 
plitude is limited to Ociit . 



with many other modes significantly affecting the fluid flow 
leading to a catastrophic decay. 

Our numerical simulations are limited to a spatial res- 
olution of Aa; ~ 0.3 km. Fluid modes with wavelengths 
smaller than A ~ 1 km cannot be accurately represented. 
It is possible that the couplings of the r-mode to a large 
number of shorter wavelength modes can limit Ocr it to a n 
even smaller valu e (as suggested by lArras et al. i2003l) : 
iBrink et al.l i2005fl 'l. The amplitude for the catastrophic de- 
cay Ocrit determined in this paper is as an upper bound. 

The study in this paper is carried out with the neutron 
star described as a self-gravitating perfect fluid body with a 
polytropic equation of state. No particular dissipation mech- 
anism is assumed. This suggests that the hydrodynamic be- 
havior we see is of rather general nature. 

Finally, we comment on future prospects of numerical 
investigations of the r-mode instability. An interesting ques- 
tion that one might ask is whether numerical simulations 
could be used to test the maximal amplitudes inferred in this 
paper by starting the r-mode evolution directly at a ~ Ocrit 
instead of a ~ 0(1) as wc have done in our simulations. As 
we have shown in|Grcssrnan ct al. (20jD^, the time it takes 
for the r-mode to reach the point of catastrophic decay de- 
pends sensitively on the value of a. For a ~ Ocrit, it would 
take a time much longer than what could be modelled accu- 
rately by our finite-differencing code in order to see a signifl- 
cant growth of the daughter modes and the catastrophic de- 
cay of the r-mode. Nevertheless, it might be possible to use 
hydrodynamics code s based on spectral methods (see, e.g., 
IVillain fc Bonazzolal i2002fl ') to handle such long timescale 
problem, since spectral methods can achieve a much higher 
numerical accuracy in general for smooth fluid flow. 



5 CONCLUSIONS AND DISCUSSIONS 

In this paper, we identify a particular three-mode coupling 
between the r-mode and the pair of fluid modes respon- 
sible for the catastrophic decay of large-amplitude r-modes 
seen in our previous simulations CGressman et al.,20o3) . The 
three modes satisfy a suitable resonance condition to high 
accuracy. This also demonstrates that hydrodynamic three- 
mode coupling can play an i mportant role in th e r-mode 
evolution as first pointed out bv lArras et al .' (2003) . We have 
also extracted the eigenfunctions of the two daughter modes, 
and demonstrated that they are consistent with having az- 
imuthal mode number ra = —1. 

While the fully nonlinear simulations are carried out 
with a large amplitude r-mode that is obtained by an ar- 
tificially large gravitational radiation reaction, we succeed 
in drawing conclusions on the r-mode evolution at smaller 
amplitude through a comparison with a second-order per- 
turbation analysis. In particular, we show that the three- 
mode equations J^-Q are valid in the simulation when 
the daughter modes are small, and determine the coupling 
strength \hb.ab\ of this particular three-mode coupling to 
be (at least) ~ O(10~*). This implies that the two daugh- 
ter modes will catch up with the r-mode at a rather small 
amplitude (a — Ocrit ~ O(10~^)), for an r-mode driven 
by gravitational radiation reaction. We conclude that at 
a ~ O(10~^) the usual picture of the fiuid flow of the neu- 
tron star dominated by an r-mode would no longer be valid. 
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APPENDIX A: TRANSFORMATION RULES OF 
TOROIDAL MODES 

In this appendix, we derive the transformation rules for 
the Cartesian velocity components of toroidal modes on the 
equatorial plane. An arbitrary fluid displacement vector ^ 
can be decomposed in vector spherical harmonics as 

i = {AlmYlmf + BlmVYlm + C(^r X VF;^) , (Al) 

Im 

where {f,9,(:f>) is the coordinate vector basis associated 
with the spherical polar coordinates {r,9,4i); Yim{d,(l)) are 
the standard spherical harmonics; the quantities Aim{r), 
Bim{r), and C;,n(') are functions of r only. Fluid modes 
for which Cim ~ (for all / and m) are called spheroidal 
modes, while those for which Aim = = are called 
toroidal modes. Here we only consider the latter class of 
modes. 

The spherical harmonics can be written as Yim = 
bim.e''"^'^ Pi^ {cos 9) , where bim is a constant depending on / 
and m, and are the associated Legendre polynomials. 
For a fixed (I, m) component, we have the following angular 
dependence for torodial modes 



-Pr{cose)9 + sin 6iPr'(cos 61)0 



(A2) 



where Pl"'(x) = dP["/dx. Considering the physically rele- 
vant case of real ^, the fluid motion of a torodial mode on 
the equatorial plane {9 = n/2) is 

6m ~ msinm<?iP('"(0)^ — cos m^Pl^' {0)4>. (A3) 

Using the standard transformation between the spherical ba- 
sis vectors (f, 9, (p) and the Cartesian one {x, y, z), wc obtain 
the following angular dependence of the Cartesian compo- 
nents of 6^ on the equatorial plane: 

~ cos m<p sin 0Pr ' (0) , ( A4) 

~ - cos mct> cos </>Pr' (0) , ( A5) 

~ -m sin m0Pr (0). (A6) 

We see that under a 7r-rotation {(j> ^ (j> + n), the Cartesian 
components transform according to 

-> (-l)'"+'6™, (A7) 

CL - (-i)'"+'CL. (A8) 

- i-irem- (A9) 



